#!/bin/sh

# run from directory where this script is
cd `dirname $0`
EXAMPLE_DIR=`pwd`

# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi

$ECHO
$ECHO "$EXAMPLE_DIR : starting"
$ECHO
$ECHO "This example shows how to calculate the phonon-induced renormalization"
$ECHO "of electron bands for diamond and how to use the Fourier interpolation"
$ECHO "of phonon potential."

# set the needed environment variables
. ../../../environment_variables

# required executables and pseudopotentials
BIN_LIST="pw.x ph.x q2r.x matdyn.x dvscf_q2r.x postahc.x"
PSEUDO_LIST="C.UPF"

$ECHO
$ECHO "  executables directory: $BIN_DIR"
$ECHO "  pseudo directory:      $PSEUDO_DIR"
$ECHO "  temporary directory:   $TMP_DIR"
$ECHO "  checking that needed directories and files exist...\c"

# check for gnuplot
GP_COMMAND=`which gnuplot 2>/dev/null`
if [ "$GP_COMMAND" = "" ]; then
        $ECHO
        $ECHO "gnuplot not in PATH"
        $ECHO "Results will not be plotted"
fi

# check for directories
for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
    if test ! -d $DIR ; then
        $ECHO
        $ECHO "ERROR: $DIR not existent or not a directory"
        $ECHO "Aborting"
        exit 1
    fi
done
for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
    if test ! -d $DIR ; then
        mkdir $DIR
    fi
done
cd $EXAMPLE_DIR/results

# check for executables
for FILE in $BIN_LIST ; do
    if test ! -x $BIN_DIR/$FILE ; then
        $ECHO
        $ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
        $ECHO "Aborting"
        exit 1
    fi
done

# check for pseudopotentials
for FILE in $PSEUDO_LIST ; do
    if test ! -r $PSEUDO_DIR/$FILE ; then
       $ECHO
       $ECHO "Downloading $FILE to $PSEUDO_DIR...\c"
            $WGET $PSEUDO_DIR/$FILE $NETWORK_PSEUDO/$FILE 2> /dev/null
    fi
    if test $? != 0; then
        $ECHO
        $ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
        $ECHO "Aborting"
        exit 1
    fi
done
$ECHO " done"

# how to run executables
PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
PH_COMMAND="$PARA_PREFIX $BIN_DIR/ph.x $PARA_POSTFIX"
MATDYN_COMMAND="$PARA_PREFIX $BIN_DIR/matdyn.x $PARA_POSTFIX"
Q2R_COMMAND="$PARA_PREFIX $BIN_DIR/q2r.x $PARA_POSTFIX"
DVSCFQ2R_COMMAND="$PARA_PREFIX $BIN_DIR/dvscf_q2r.x $PARA_POSTFIX"
POSTAHC_COMMAND="$PARA_PREFIX $BIN_DIR/postahc.x $PARA_POSTFIX"
$ECHO
$ECHO "  running pw.x as:     $PW_COMMAND"
$ECHO "  running ph.x as:     $PH_COMMAND"
$ECHO "  running q2r.x as:    $Q2R_COMMAND"
$ECHO "  running matdyn.x as: $MATDYN_COMMAND"
$ECHO "  running dvscf_q2r.x as: $DVSCFQ2R_COMMAND"
$ECHO "  running postahc.x as: $POSTAHC_COMMAND"
$ECHO

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/diam*
rm -rf $TMP_DIR/_ph0/diam*
$ECHO " done"

# self-consistent calculation
cat > diam.scf.in << EOF
 &control
    calculation     = 'scf'
    prefix          = 'diam'
    restart_mode    = 'from_scratch'
    pseudo_dir      = '$PSEUDO_DIR/'
    outdir          = '$TMP_DIR/'
    tprnfor         = .true.
    tstress         = .true.
 /
 &system
    ibrav           = 2
    celldm(1)       = 6.64
    nat             = 2
    ntyp            = 1
    ecutwfc         = 60
    occupations     = 'fixed'
 /
 &electrons
    diagonalization = 'david'
    mixing_beta     = 0.7
    conv_thr        = 1.0d-12
 /
ATOMIC_SPECIES
  C  12.01078  C.UPF
ATOMIC_POSITIONS alat
  C   0.00  0.00  0.00
  C   0.25  0.25  0.25
K_POINTS automatic
 6 6 6 0 0 0
EOF
$ECHO "  running the scf calculation...\c"
$PW_COMMAND < diam.scf.in > diam.scf.out
check_failure $?
$ECHO " done"

# phonon calculation on a (444) uniform grid of q-points
cat > diam.ph.in << EOF
phonons of diamond on a coarse 333 grid
&inputph
  prefix   = 'diam'
  outdir   = '$TMP_DIR/'
  reduce_io = .true.
  fildyn   = 'diam.dyn'
  fildvscf = 'dvscf'
  ldisp    = .true.
  nq1 = 3
  nq2 = 3
  nq3 = 3
  tr2_ph   =  1.0d-16
 /
EOF
$ECHO "  running the phonon calculation ...\c"
$PH_COMMAND < diam.ph.in > diam.ph.out
check_failure $?
$ECHO " done"

# q2r
cat > q2r.in <<EOF
&input
  fildyn = 'diam.dyn'
  zasr   = 'crystal'
  flfrc  = 'diam.ifc'
/
EOF

$ECHO "  transforming C(q) => C(R) ...\c"
$Q2R_COMMAND < q2r.in > q2r.out
check_failure $?
$ECHO " done"

# dvscf_q2r
cat > dvscfq2r.in <<EOF
&input
  prefix = 'diam'
  outdir = '$TMP_DIR/'
  fildyn = 'diam.dyn'
  fildvscf = 'dvscf'
  wpot_dir = 'wpot/'

  do_long_range = .false.
  do_charge_neutral = .false.
 /
EOF

$ECHO "  transforming phonon potential V(q) => V(R) ...\c"
$DVSCFQ2R_COMMAND < dvscfq2r.in > dvscfq2r.out
check_failure $?
$ECHO " done"

# run SCF again
$ECHO "  running the scf calculation again ...\c"
$PW_COMMAND < diam.scf.in > diam.scf.out
check_failure $?
$ECHO " done"

# NSCF
cat > diam.nscf.in << EOF
 &control
    calculation     = 'nscf'
    prefix          = 'diam'
    restart_mode    = 'from_scratch'
    pseudo_dir      = '$PSEUDO_DIR/'
    outdir          = '$TMP_DIR/'
 /
 &system
    ibrav           = 2
    celldm(1)       = 6.64
    nat             = 2
    ntyp            = 1
    ecutwfc         = 60
    occupations     = 'fixed'
    nbnd = 15
    nosym = .true.
    noinv = .true.
 /
 &electrons
    diago_full_acc  = .true.
    diagonalization = 'david'
    mixing_beta     = 0.7
    conv_thr        = 1.0d-10
 /
ATOMIC_SPECIES
  C  12.01078  C.UPF
ATOMIC_POSITIONS alat
  C   0.00  0.00  0.00
  C   0.25  0.25  0.25
K_POINTS crystal
     2
  0.000  0.000  0.000    1.0
  0.365  0.365  0.000    1.0
EOF
$ECHO "  running the nscf calculation ...\c"
$PW_COMMAND < diam.nscf.in > diam.nscf.out
check_failure $?
$ECHO " done"

# ph.x ahc with coarse grid
cat > diam.ahc.coarse.in << EOF
electron-phonon of diamond on a coarse 333 grid
&inputph
  prefix   = 'diam'
  outdir   = '$TMP_DIR/'
  reduce_io = .true.
  fildyn   = 'dyn_dir_ahc_coarse/diam.dyn'
  tr2_ph   =  1.0d-20

  ldisp    = .true.
  nq1      = 3
  nq2      = 3
  nq3      = 3

  ldvscf_interpolate = .true.
  wpot_dir = 'wpot/'

  trans    = .false.
  electron_phonon = 'ahc'
  ahc_nbnd = 8

  ahc_dir = 'ahc_dir_coarse'
 /
EOF
mkdir dyn_dir_ahc_coarse
$ECHO "  running the e-ph calculation on a coarse (333) q-grid ...\c"
$PH_COMMAND < diam.ahc.coarse.in > diam.ahc.coarse.out
check_failure $?
$ECHO " done"

# matdyn for coarse grid
cat > matdyn.coarse.in <<EOF
 &input
    asr = 'crystal'
    amass(1) = 12.01078
    flfrc = 'diam.ifc'
    flvec = 'diam.modes_coarse'
    flfrq = ''
    fleig = ''
    q_in_band_form = .false.
    q_in_cryst_coord = .false.

    loto_disable = .true.
 /
  27
   0.000000000000000E+00   0.000000000000000E+00   0.000000000000000E+00
  -0.333333333333333E+00   0.333333333333333E+00  -0.333333333333333E+00
   0.333333333333333E+00  -0.333333333333333E+00   0.333333333333333E+00
   0.333333333333333E+00   0.333333333333333E+00   0.333333333333333E+00
   0.000000000000000E+00   0.666666666666667E+00   0.000000000000000E+00
   0.666666666666667E+00  -0.555111512312578E-16   0.666666666666667E+00
  -0.333333333333333E+00  -0.333333333333333E+00  -0.333333333333333E+00
  -0.666666666666667E+00  -0.555111512312578E-16  -0.666666666666667E+00
   0.000000000000000E+00  -0.666666666666667E+00   0.000000000000000E+00
  -0.333333333333333E+00  -0.333333333333333E+00   0.333333333333333E+00
  -0.666666666666667E+00   0.000000000000000E+00   0.000000000000000E+00
   0.555111512312578E-16  -0.666666666666667E+00   0.666666666666667E+00
   0.000000000000000E+00   0.000000000000000E+00   0.666666666666667E+00
  -0.333333333333333E+00   0.333333333333333E+00   0.333333333333333E+00
   0.333333333333333E+00  -0.333333333333333E+00   0.100000000000000E+01
  -0.666666666666667E+00  -0.666666666666667E+00  -0.555111512312578E-16
  -0.100000000000000E+01  -0.333333333333333E+00  -0.333333333333333E+00
  -0.333333333333333E+00  -0.100000000000000E+01   0.333333333333333E+00
   0.333333333333333E+00   0.333333333333333E+00  -0.333333333333333E+00
   0.555111512312578E-16   0.666666666666667E+00  -0.666666666666667E+00
   0.666666666666667E+00   0.000000000000000E+00   0.000000000000000E+00
   0.666666666666667E+00   0.666666666666667E+00  -0.555111512312578E-16
   0.333333333333333E+00   0.100000000000000E+01  -0.333333333333333E+00
   0.100000000000000E+01   0.333333333333333E+00   0.333333333333333E+00
   0.000000000000000E+00   0.000000000000000E+00  -0.666666666666667E+00
  -0.333333333333333E+00   0.333333333333333E+00  -0.100000000000000E+01
   0.333333333333333E+00  -0.333333333333333E+00  -0.333333333333333E+00
EOF

$ECHO "  calculating omega(q) from C(R) for a coarse q-grid ...\c"
$MATDYN_COMMAND < matdyn.coarse.in > matdyn.coarse.out
check_failure $?
$ECHO " done"

# postahc.x with coarse grid
cat > postahc.coarse.in << EOF
&input
  prefix = 'diam'
  outdir = '$TMP_DIR/'

  ahc_dir = 'ahc_dir_coarse/'
  flvec = 'diam.modes_coarse'
  ahc_nbnd = 8
  amass_amu(1) = 12.01078
  ahc_nbndskip = 0

  eta_eV = 0.1
  efermi_eV = 16.0
  temp_kelvin = 300.0
/
EOF
$ECHO "  calculating electron self-energy with a coarse q-grid ...\c"
$POSTAHC_COMMAND < postahc.coarse.in > postahc.coarse.out
mv selfen_real.dat selfen_coarse_real.dat
mv selfen_imag.dat selfen_coarse_imag.dat
check_failure $?
$ECHO " done"

# ph.x ahc with fine grid
cat > diam.ahc.fine.in << EOF
electron-phonon of diamond on a fine 444 grid
&inputph
  prefix   = 'diam'
  outdir   = '$TMP_DIR/'
  reduce_io = .true.
  fildyn   = 'dyn_dir_ahc_fine/diam.dyn'
  tr2_ph   =  1.0d-20

  ldisp    = .true.
  nq1      = 4
  nq2      = 4
  nq3      = 4

  ldvscf_interpolate = .true.
  wpot_dir = 'wpot/'

  trans    = .false.
  electron_phonon = 'ahc'
  ahc_nbnd = 8

  ahc_dir = 'ahc_dir_fine'
  skip_upper = .true.
 /
EOF
mkdir dyn_dir_ahc_fine
$ECHO "  running the e-ph calculation on a fine (444) q-grid ...\c"
$PH_COMMAND < diam.ahc.fine.in > diam.ahc.fine.out
check_failure $?
$ECHO " done"

# matdyn for fine grid
cat > matdyn.fine.in <<EOF
 &input
    asr = 'crystal'
    amass(1) = 12.01078
    flfrc = 'diam.ifc'
    flvec = 'diam.modes_fine'
    flfrq = ''
    fleig = ''
    q_in_band_form = .false.
    q_in_cryst_coord = .false.

    loto_disable = .true.
 /
  64
   0.000000000000000E+00   0.000000000000000E+00   0.000000000000000E+00
  -0.250000000000000E+00   0.250000000000000E+00  -0.250000000000000E+00
   0.500000000000000E+00  -0.500000000000000E+00   0.500000000000000E+00
   0.250000000000000E+00  -0.250000000000000E+00   0.250000000000000E+00
   0.250000000000000E+00   0.250000000000000E+00   0.250000000000000E+00
   0.000000000000000E+00   0.500000000000000E+00   0.000000000000000E+00
   0.750000000000000E+00  -0.250000000000000E+00   0.750000000000000E+00
   0.500000000000000E+00   0.000000000000000E+00   0.500000000000000E+00
  -0.500000000000000E+00  -0.500000000000000E+00  -0.500000000000000E+00
  -0.750000000000000E+00  -0.250000000000000E+00  -0.750000000000000E+00
   0.000000000000000E+00  -0.100000000000000E+01   0.000000000000000E+00
  -0.250000000000000E+00  -0.750000000000000E+00  -0.250000000000000E+00
  -0.250000000000000E+00  -0.250000000000000E+00  -0.250000000000000E+00
  -0.500000000000000E+00   0.000000000000000E+00  -0.500000000000000E+00
   0.250000000000000E+00  -0.750000000000000E+00   0.250000000000000E+00
   0.000000000000000E+00  -0.500000000000000E+00   0.000000000000000E+00
  -0.250000000000000E+00  -0.250000000000000E+00   0.250000000000000E+00
  -0.500000000000000E+00   0.000000000000000E+00   0.000000000000000E+00
   0.250000000000000E+00  -0.750000000000000E+00   0.750000000000000E+00
   0.000000000000000E+00  -0.500000000000000E+00   0.500000000000000E+00
   0.000000000000000E+00   0.000000000000000E+00   0.500000000000000E+00
  -0.250000000000000E+00   0.250000000000000E+00   0.250000000000000E+00
   0.500000000000000E+00  -0.500000000000000E+00   0.100000000000000E+01
   0.250000000000000E+00  -0.250000000000000E+00   0.750000000000000E+00
  -0.750000000000000E+00  -0.750000000000000E+00  -0.250000000000000E+00
  -0.100000000000000E+01  -0.500000000000000E+00  -0.500000000000000E+00
  -0.250000000000000E+00  -0.125000000000000E+01   0.250000000000000E+00
  -0.500000000000000E+00  -0.100000000000000E+01   0.000000000000000E+00
  -0.500000000000000E+00  -0.500000000000000E+00   0.000000000000000E+00
  -0.750000000000000E+00  -0.250000000000000E+00  -0.250000000000000E+00
   0.000000000000000E+00  -0.100000000000000E+01   0.500000000000000E+00
  -0.250000000000000E+00  -0.750000000000000E+00   0.250000000000000E+00
   0.500000000000000E+00   0.500000000000000E+00  -0.500000000000000E+00
   0.250000000000000E+00   0.750000000000000E+00  -0.750000000000000E+00
   0.100000000000000E+01   0.000000000000000E+00   0.000000000000000E+00
   0.750000000000000E+00   0.250000000000000E+00  -0.250000000000000E+00
   0.750000000000000E+00   0.750000000000000E+00  -0.250000000000000E+00
   0.500000000000000E+00   0.100000000000000E+01  -0.500000000000000E+00
   0.125000000000000E+01   0.250000000000000E+00   0.250000000000000E+00
   0.100000000000000E+01   0.500000000000000E+00   0.000000000000000E+00
   0.000000000000000E+00   0.000000000000000E+00  -0.100000000000000E+01
  -0.250000000000000E+00   0.250000000000000E+00  -0.125000000000000E+01
   0.500000000000000E+00  -0.500000000000000E+00  -0.500000000000000E+00
   0.250000000000000E+00  -0.250000000000000E+00  -0.750000000000000E+00
   0.250000000000000E+00   0.250000000000000E+00  -0.750000000000000E+00
   0.000000000000000E+00   0.500000000000000E+00  -0.100000000000000E+01
   0.750000000000000E+00  -0.250000000000000E+00  -0.250000000000000E+00
   0.500000000000000E+00   0.000000000000000E+00  -0.500000000000000E+00
   0.250000000000000E+00   0.250000000000000E+00  -0.250000000000000E+00
   0.000000000000000E+00   0.500000000000000E+00  -0.500000000000000E+00
   0.750000000000000E+00  -0.250000000000000E+00   0.250000000000000E+00
   0.500000000000000E+00   0.000000000000000E+00   0.000000000000000E+00
   0.500000000000000E+00   0.500000000000000E+00   0.000000000000000E+00
   0.250000000000000E+00   0.750000000000000E+00  -0.250000000000000E+00
   0.100000000000000E+01   0.000000000000000E+00   0.500000000000000E+00
   0.750000000000000E+00   0.250000000000000E+00   0.250000000000000E+00
  -0.250000000000000E+00  -0.250000000000000E+00  -0.750000000000000E+00
  -0.500000000000000E+00   0.000000000000000E+00  -0.100000000000000E+01
   0.250000000000000E+00  -0.750000000000000E+00  -0.250000000000000E+00
   0.000000000000000E+00  -0.500000000000000E+00  -0.500000000000000E+00
   0.000000000000000E+00   0.000000000000000E+00  -0.500000000000000E+00
  -0.250000000000000E+00   0.250000000000000E+00  -0.750000000000000E+00
   0.500000000000000E+00  -0.500000000000000E+00   0.000000000000000E+00
   0.250000000000000E+00  -0.250000000000000E+00  -0.250000000000000E+00
EOF

$ECHO "  calculating omega(q) from C(R) for a fine q-grid ...\c"
$MATDYN_COMMAND < matdyn.fine.in > matdyn.fine.out
check_failure $?
$ECHO " done"

# postahc.x with fine grid
cat > postahc.fine.in << EOF
&input
  prefix = 'diam'
  outdir = '$TMP_DIR/'

  ahc_dir = 'ahc_dir_fine/'
  flvec = 'diam.modes_fine'
  ahc_nbnd = 8
  ahc_nbndskip = 0

  eta_eV = 0.1
  efermi_eV = 16.0
  temp_kelvin = 300.0

  skip_upper = .true.
  skip_dw = .true.

  amass_amu(1) = 12.01078
/
EOF
$ECHO "  calculating electron self-energy with a fine q-grid ...\c"
$POSTAHC_COMMAND < postahc.fine.in > postahc.fine.out
mv selfen_real.dat selfen_fine_real.dat
mv selfen_imag.dat selfen_fine_imag.dat
check_failure $?
$ECHO " done"

$ECHO
$ECHO "$EXAMPLE_DIR: done"
